Notice that the Euler method is actually equivalent to a sequence of isogeometric L2 projections, similar to the calculation of a bitmap projection. This time, however, our "bitmap" is the state of the system in the previous time step, plus changes caused by the "physics" of the modeled phenomenon during the time step, plus changes caused by the force acting on our system during the time step.
The method of treating the non-stationary equation as sequences of isogeometric projections is described in [1].
To develop a solver using the Euler method in the isogeometric finite element method, we need to transform a strong formulation into a weak formulation.
So we multiply our equation by the test functions
\( (u_{t+1},w) = (u_t + dt * \mathcal{L}(u_t)+ dt* f_t,w) \)
We will use a linear combination of the B-spline function to approximate the state of our system at a given moment in time. For this purpose, we select the basis of two-dimensional B-spline functions, specifying the node vectors in the direction of the axis of the coordinate system. To establish attention, we can choose a two-dimensional basis of the second order B-spline function
\( \{B^x_{i,2}(x)B^y_{j,2}(y)\}_{i=1,...,N_x;j=1,...,N_y} \).
They will be used to approximate the simulated scalar field at given time moment
\( u(x,y;t+1)\approx \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j } B^x_i(x)B^y_j(y) \).
Similarly, for testing, we will use the B-spline basis functions:
\( \{ B^x_k(x)B^y_l(y) \}_{k=1,...,N_x;l=1,...,N_y } \)
Our equation now looks identical to the isgeometric L2 projection problem
\( \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j } (B^x_i(x)B^y_j(y),B^x_k(x)B^y_l(y))=RHS \quad \forall k=1,...,N_x; l=1,...,N_y \)
Our right side is not a bitmap, however, but a projection of the sum of three elements:

1. The state of our system at the previous moment in time \( (u_t,w)= \sum_{i=1,...,N_x;j=1,...,N_y} u^t_{i,j } (B^x_i(x)B^y_j(y),B^x_k(x)B^y_l(y)) \) (also multiplied by the test function and integrated over the area). Of course, we also use a linear combination of B-spline basis functions to represent the state of our system in the previous time step \( u(x,y;t)=u_t=\sum_{i=1,...,N_x;j=1,...,N_y} u^t_{i,j } B^x_i(x)B^y_j(y) \)
2. Changes induced by the "physics" of the simulated phenomenon during the time step. These changes are calculated by applying the differential operator representing the modeled phenomenon to the state of the system at the previous moment. The state of the system is of course represented by a linear combination of the B-spline function. So the differential operator is applied to the B-spline function \( (dt * \mathcal{L}(u_t),w)=\sum_{i=1,...,N_x;j=1,...,N_y} dt*u^{t}_{i,j } ({\cal L} (B^x_i(x)B^y_j(y)),B^x_k(x)B^y_l(y)) \). For the heat transfer problem, we have \( \sum_{i=1,...,N_x;j=1,...,N_y} dt*u^{t}_{i,j } (\frac{\partial^2 (B^x_i(x)B^y_j(y))}{\partial x^2}+\frac{\partial^2 (B^x_i(x)B^y_j(y))}{\partial y^2},B^x_k(x)B^y_l(y)) \). We can integrate by parts \( \sum_{i=1,...,N_x;j=1,...,N_y} dt*u^{t}_{i,j } \left( (\frac{\partial (B^x_i(x)B^y_j(y))}{\partial x}B^y_j(y),\frac{\partial (B^x_k(x)B^y_l(y))}{\partial x}B^y_l(y))+(B^x_i(x)\frac{\partial (B^x_i(x)B^y_j(y))}{\partial y},B^x_k(x)\frac{\partial (B^x_k(x)B^y_l(y))}{\partial y})\right) \) which due to the structure of the tensor product of the B-spline function, and due to the fact that the derivative in the y direction of the B-spline function of the variable x gives 0, and the derivative in the x direction of the B-spline function of the variable y gives 0, we get \( \sum_{i=1,...,N_x;j=1,...,N_y} dt*u^{t}_{i,j } \left( (\frac{\partial B^x_i(x)}{\partial x}B^y_j(y),\frac{\partial B^x_k(x)}{\partial x}B^y_l(y))+(B^x_i(x)\frac{\partial B^y_j(y)}{\partial y},B^x_k(x)\frac{\partial B^y_l(y)}{\partial y})\right) \).
3. Changes caused by the force acting on the system during the time step \( (f_t,w)= (f_t(x,y),B^x_k(x)B^y_l(y)) \)

In the explicit method, the system matrix is a mass matrix \( \{ M_x \otimes M_y \}_{i,j,k,l}= \int B^x_i(x)B^x_k(x)B^y_j(y)B^y_l(y) = \int B^x_i(x)B^y_j(y)B^x_k(x)B^y_l(y) = {\bf M}_{i,j,k,l } \)
which has the structure of a Kronecker product and can be factored in linear time by an alternating-direction solver. The disadvantage of the method is explicitly the fact that a single time step cannot be too long, otherwise the method will be unstable. It will show up with an explosion of numerical simulations. This is due to the so-called Courant-Friedrichs-Lewy condition (CFL condition). In general, it says numerical simulation in which we model a scalar field \( u(x,t) \) is stable if \( C= \frac{u * dt}{ h} < C_{max } \) where \( u = max u(x,t) \) denotes the maximum value of the simulated field, \( dt \) denotes the time step, \( h \) denotes element dimension (which in case of IGA translates to the knot points span), and \( C_{max} \) is a problem specific constant. The practical significance of this condition is as follows: if we increase the number of mesh elements (increase the number of points in the knot vector) then \( h \) will decrease. To make sure our simulation does not "explode," we have to decrease the time step. This is a big limitation because it means that increasing accuracy in space forces time step sizes to decrease, which means a higher simulation cost. One solution to this problem is the use of implicit simulations with time integration schemes that allow the use of alternating direction solvers.